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Abstract. The dynamical, long-wavelength longitudinal and transverse exchange- 
correlation potentials for a homogeneous electron gas are evaluated in a microscopic 
model based on an approximate decoupling of the equation of motion for the current- 
current response function. The transverse spectrum turns out to be very similar to the 
longitudinal one. We obtain evidence for a strong spectral structure near twice the 
^ ' plasma frequency due to a two-plasmon threshold for two— pair excitations, which 

(N ■ may be observable in inelastic scattering experiments. Our results give the entire 

, input needed to implement the Time-Dependent Current Density Functional Theory 

C " 3 , scheme recently developed by G. Vignale and W. Kohn [Phys. Rev. Lett. 77, 2037 

1 (1996)] and are fitted to analytic functions to facilitate such applications. 

o 
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PACS numbers: 71.45Gm, 36.40+d, 78.30.Fs, 73.20.Mf 
J2 I The Time-Dependent Density Functional Theory (TD-DFT) of Runge and Gross 

101 in principle allows the study of dynamical properties of interacting many- 
^ , particle systems, which are not addressed by static DFT. The adiabatic local 

■ density approximation (ALDA) Q has enabled successful application of the theory 
Q , in the low-frequency limit. The search for a local dynamical approximation to the 
O ■ exchange-correlation (xc) potential has been frustrated by the appearence of various 

inconsistencies [H) , which can be tracked down to the non-existence of a gradient 
expansion for the frequency-dependent xc potential in terms of the density alone. 
. , Recently Vignale and Kohn (VK) have shown how a local approximation for the xc 

■ vector potential a^'^(r, w) can be obtained within linearized Current Density Functional 
Theory |^. They prove that for slowly varying densities known symmetries 
and conservation laws allow one to express a^'^(r, w) exactly in terms of the long- 
wavelength xc potentials f^c'^iuj) for a homogeneous electron gas, 

ar(r,c^) = - ^ {a.[/i;V • (nov) - 5f^,v ■ Vno] + 5f^,id,no)V ■ v 

+/Jc h»^o(V X V X v), + 2{djVi + diWj)djno - 4(a,no)V • v] 
+ no (9»v, + 9,v,) - 2(9,/JjV • v) } . (1) 

Here, no(r) is the unperturbed ground-state density, v(r,ti;) — ji(r,ti;)/no(r) is the 
local velocity and (5/^^(a;,n) = f^^{uj,n) — /^^(0,n). The potentials Z^^'^, which are 
functions of uj and of the local density no{r), are defined as the fc — ^ limit of 

f^/{k, C.) = I \^ , , 1 - VL.T (2) 

^ Xlt(^'^) + /o/™ XL,T{k,uj) + p/m \ 
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where Xl(t) is the current-current longitudinal (transverse) response function of the 
homogeneous system at density p = rio(r), Xi,(t) ^^"^ equivalent ideal-gas quantity, 
vl = 47re^/A:^ and vt = 0. 

The quest for a frequency-dependent local field factor G{k,uo) — ~ f^^{k,uo) / vl 
for the density-density response is a long-lasting problem in electron gas theory 

^, |l^, pi], p2| , p^ , p^ , [T^ . Its asymptotic behaviours are known from sum rule 
arguments [^,^ and from second-order perturbative expansions [O, ll^, and 
a smooth interpolation scheme was proposed by Gross and Kohn |H nQ] . A fully 
microscopic result, derived from an approximate treatment of two-pair excitations, 
was obtained by Bohm et al p6| , p7| using an expression first derived by Hasegawa 
and Watabe The longitudinal component of the current-current xc potential 

f^Juj) can be obtained from such results, but no sensible approximations have been 
proposed for the transverse component, preventing implementation of the Vignale- 
Kohn theory. 

In this letter we develop a full treatment of the two-pair excitation spectrum, 
including its transverse component. Our results give the entire input necessary to use 
the VK expression. 

The xc potential can be obtained within the equation of motion formalism. A 
general response function is defined as 



« A; B »^ = -I / dte'('^+")* < [A{t),B{0)] > 



(3) 



where A{t) = e'^*^e~*^* and < . 
current-current response Xij{k,^) 



> denotes ground state expectation values. The 
satisfies the equation of motion 



(4) 



The last term yields four-point response functions ^ AB; CD involving pj^ and 
jj^ operators and can be approximately decoupled by an RPA-like scheme: 



Im < AB: CD 



[Im< A;C>, 

Im< A; >, 



Im< > 
Im< B;C> 

]■ (5) 



Such a decoupling includes by construction two-pair processes, which are the lowest 
order processes with non-zero spectral strength in the relevant region of the 
plane. Perturbative approaches evaluate the LHS of equation (5) for an ideal gas and 
obtain a spectrum restricted to single pair excitations, which is essentially equivalent 
to the one obtained by using non-interacting response functions in the RHS (the 
difference lies in exchange processes, as discussed below). We intend to include the 
effect of plasmons and therefore use RPA response functions in the RHS. 

Quite lengthy calculations lead in this approach to the following result for the xc 
potentials f^d'^i^) terms of the response functions introduced in equation (0): 



d^q 



TT J (27r)3p2 1{oJ-Uj')^ 



ImxL(g,t^-w') 



aL,T—p^lmxL{q,uj') 



bL,T^lrnxTiq,^') 



(6) 



with cl ~ 23/30, ax ~ 8/15, = 8/15 and — 2/5. The expression for the 
longitudinal part is equivalent to the one obtained by Hasegawa and Watabe |ll| by 
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Table 1. Exact limiting behaviours of fxi'^i'^) from the Monte Carlo equation of 
state, in units of 2u!pi/p, and best fit parameters for Im/^^. according to (12). 





fLio) 




fLioo) 


/3 


lO^co 


lO^ci 




UJ2 


do 


lO^di 


0.5 


-0.04246 


-0.01794 


0.0177 


1.87 


0.175 


0.694 


1.75 


-3.59 


0.173 


5.72 


1 


-0.0611 


-0.0216 


0.0284 


1.48 


0.421 


1.76 


0.982 


-1.45 


0.291 


9.38 


2 


-0.0891 


-0.0252 


0.0457 


1.22 


0.895 


3.87 


0.347 


0.181 


0.49 


13.2 


3 


-0.1119 


-0.0280 


0.0600 


1.1 


1.29 


6.09 


0.143 


0.693 


0.664 


16.7 


4 


-0.1320 


-0.0308 


0.0724 


1.02 


1.65 


7.87 


-0.143 


1.33 


0.824 


17 


5 


-0.1503 


-0.0338 


0.0835 


0.955 


1.94 


9.82 


-0.27 


1.61 


0.974 


18.3 


6 


-0.1674 


-0.0370 


0.0935 


0.899 


2.22 


11.6 


-0.361 


1.82 


1.12 


19.3 


10 


-0.2276 


-0.0518 


0.1267 


0.698 


3.11 


17.9 


-0.565 


2.27 


1.64 


22.1 


15 


-0.2917 


-0.0725 


0.1587 


0.474 


3.94 


24.2 


-0.69 


2.54 


2.22 


23.9 


20 


-0.3483 


-0.0939 


0.1847 


0.259 


5.54 


24.7 


-0.808 


2.78 


2.75 


22.8 



diagrammatic means, and similar to the one obtained by Neilson et al. |15|| (see the 
discussion in |l6[|). The resuh for the transverse component is new. 

The real part of fxc can be obtained via the Kramers-Kronig relations, namely 

Re/,,(.;) ^ Re/,,(oo) + - d^'I^^^M^^ . (7) 



fixes the high-frequency limit 

[dL.T{{ke)^{ke)°) + eL,T{pe)] 



(8) 



The first term in the RHS of equation 

1 

Tp 

with (11 = 4, dr = 4/3, cl ~ 8/15 and er = —4/15. The average kinetic and potential 
energies (fee) and (pe) can be obtained from the Monte Carlo equation of state jist , 
as can the compressibility Kt which fixes the longitudinal static limit 
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(9) 



In the above equations, (fce)° and Kj. denote ideal gas results. The asymptotic values 
of /Jc"^ are listed in Table |l} We remark that there is no proof that the order of 
limits in equation (^) can be safely interchanged nor in general that f^^{k,uj) is a 
continuous function in the limit {k^uo) — > (0,0). Indeed, there are arguments which 
indicate a small discontinuity JlOf . In view of (i) the small estimated discontinuity, 
(ii) the uncertainty in its precise value and (iii) the appeal of a theory which reduces 
continuously to LDA as a; — > 0, we prefer to enforce continuity as explained below. 

The above approximation neglects exchange processes, which in perturbative 
treatments reduce the total two-pair spectral weight by a factor of 2 if w 3> Ef, 
£f being the Fermi energy. We approximately include exchange by multiplying Imf^^ 
by the phenomenological factor 

(3 + 0.5uj/2eF 



gx{uj) 



(10) 



l + uj/2eF ' 

where the parameter (3 is used to enforce continuity aX k — oj = and turns out to be 
close to 1 at metallic densities (see Table |l]) . For consistency we use the same factor 
gx{ijj) in the transverse component. Our final high-frequency result is 

3/2 



,Ry 



(11) 
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where cl — 237r/15 and ct = 167r/15, in agreement with the result of Ghck and Long 
p3[ for the fongitudinal term and its extension to the transverse one (unpubHshed) . 

Equations (6) could be considered as the k = component of a self-consistency 
condition, but the self-consistent problem cannot be solved since the equations 
determining the finite-A: behaviour of fxc are unknown. We have evaluated equation 
(6) by using the RPA response functions in the LHS. The imaginary part of the 
RPA longitudinal susceptibility consists of two separate contributions: the first arises 
from the broad single-pair continuum, the other from the sharply peaked plasmon 
excitation. In turn this leads to a structure in \m.fxc{^) around twice the plasma 
frequency (for details see Ref. |l6|). As we are neglecting retardation the transverse 
susceptibility contains only a broad continuum and therefore has larger spectral 
strength than the longitudinal one at low frequency. 

Figures |^ and |^ report our results for the imaginary and real parts of f^c and 
compare them with the interpolation scheme of Gross and Kohn Q for the longitudinal 
term. Both curves reproduce the asymptotic limits (8) and as well as the 
^-3/2 high-frequency behaviour, but the behaviours at intermediate frequencies are 
strikingly different. Our curves for Tie fxc exhibit a sharp minimum around 2^^;, 
which corresponds to the sharp structure found in Imfxc at the same frequency. The 
physical origin lies in the large spectral strength of the plasmon excitation as compared 
to single-pair excitations, which accumulates most of the spectral strength of two-pair 
processes near 2ujpi . This spectral structure becomes sharper with increasing coupling 
strength, i.e. with increasing rg (r^ is defined as (47rpa'Q/3)'^/'^, as being the effective 
Bohr radius). We note that within the accuracy of the present model /^^(w = 0) is 
indistinguishable from zero in the entire density range explored here. 

Figure ^ shows stronger xc effects in the real part at intermediate frequencies 
(w ~ ujpi) than in both the uj = and the lu = oo limit, at variance from the Gross- 
Kohn interpolation. Available experimental data on the plasmon dispersion coefficient 
in alkali metals |Q indicate much stronger deviations from RPA than predicted by 
existing electron gas theories. Although band-structure effects play a significant role 
in a quantitative comparison, the dynamic correction indicated by our results improves 
the agreement with the experimental data. 

Figure || displays the dynamic structure factor S{k,uj) = '-{2k^/puj'^)ImxL{k,uj) 
which is relevant to inelastic scattering experiments. The threshold behaviour at 
frequency 2ujpi is a clearcut signature of the present results for fxc- 

While the detailed shapes of our results are largely approximation-dependent, we 
do not expect the main qualitative features to be substantially modified in a more 
refined theory. Both self-consistent SPA results jl^ and preliminary calculations 
with an STLS response function ||^ show the same qualitative behaviours at values of 
Ts where the plasmon dispersion coefRcient is positive. A much deeper minimum in 
Tlefxc is found at negative plasmon dispersion, which in STLS is the case for rg > 5. 

We also remark that the shapes of the longitudinal and transverse spectra are 
very similar. Equation (6) shows that Im/J^ = 16/23 • Imf^^ for u; ^ ujpi whereas 
Im/^ = 3/4 -Im/^^ for w <C Wpj. In fact, the transverse spectrum is rather accurately 
reproduced at all frequencies by setting Im/J^(a;) ~ 0.72 • Im/J"^(a;), whereas for the 
real part there is an additional shift due to the different a; = oo value. 

Our numerical results for Im/^^(a;) can be accurately reproduced by the expression 
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Figure 1. Imaginary parts of f^^(u}) and f'^^(uj) in units of 2u}pi/p, as functions of 
Lo/ujpl at rs = 1 (left panel) and rs = 5 (right panel). The present results (full curves) 
are compared with the Gross— Kohn interpolation for the longitudinal component 
(dashed curves) and with the fit discussed in the text (dotted curves). 
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Figure 3. Dynamic structure factor S{k,ui) at = 5 as a function of uj/ujpi at 
various values of krsas on a semilogarithmic scale. 
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^0^ + ^1 ^(7^-5^1 ^^^^^ 

Im/i;H = -g.H <( ^ (12) 

—( y= 7 > 2) 

where oj is in units of uspi, f^^ in units of luj^ij p. The fit parameters given in Table 
^ were obtained by imposing (i) continuity at cj = 2a;p;, (ii) conservation of the 
normalization J luifxcduj /uj and (iii) the asymptotic behaviour given by equation (jTl]). 
The remaining three parameters were fitted to the numerical data. The real part can 
be obtained from the Kramers-Kronig relation and the high-frequency values given 
in Table |l|; the low-frequency longitudinal values from the same Table can be used for 
a check. The good quality of the resulting fits is shown in Figures |and|| 

In summary, we have presented a microscopic model for the longitudinal and 
transverse xc potentials of the 3D homogeneous electron gas. We have given numerical 
evaluations yielding results which are significantly different from previous interpolation 
schemes and compatible with experimental data on the plasmon dispersion coefficient. 
Our results have been fitted to analytic expressions and provide all the input needed 
for TD-DFT computations in the linearized long-wavelength regime. 
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